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The interplay of exchange correlations and spin-orbit interaction (SOI) on the many-body spec¬ 
trum of a copper phtalocyanine (CuPc) molecule and their signatures in transport are investigated. 

We first derive a minimal model Hamiltonian in a basis of frontier orbitals which is able to re¬ 
produce experimentally observed singlet-triplet splittings; in a second step SOI effects are included 
perturbatively. Major consequences of the SOI are the splitting of former degenerate levels and a 
magnetic anisotropy, which can be captured by an effective low-energy spin Hamiltonian. We show 
that STM-based magnetoconductance measurements can yield clear signatures of both these SOI 
induced effects. 


I. INTRODUCTION 

Spin-orbit interaction (SOI) can play a major role in 
molecular spintronics. For example, in combination with 
the configuration of the non-magnetic component (or¬ 
ganic ligand), it is known to be essential in establish¬ 
ing magnetic anisotropy in high-spin molecular magnet^. 
Effective spin-Hamiltonians are commonly used to de¬ 
scribe this anisotropy, and usually well capture the low 
energy properties of these systems, see e.g. Ref.l^. Such 
effective Hamiltonians have been derived microscopically 
for widely studied molecular magnets like Feg, Fe4 and 
Mni2P. Recently, magnetic anisotropy effects could be 
directly probed by magnetotransport spectroscopy for 
Fe4 in quantum dot setupJ^. An interesting question 
is hence if other classes of metallorganic compounds, 
like the widely studied metal phthalocyanine^^, exhibit 
magnetic anisotropy induced by the interplay of elec¬ 
tronic correlations and SOI. Indeed, in an XMCD anal¬ 
ysis copper phthalocyanine (CuPc) was found to exhibit 
enormous anisotropies in both spin and orbital dipole mo- 
mentJ^. Furthermore, recent experimental findings for 
cobalt ptht halo cyanine in an STM setupP suggest that 
many-body correlations play an important role in the in¬ 
terpretation of the transport measurements. In recent 
worlP, we have explictly investigated long range and 
short range electron-electron correlations effects in CuPc 
and found a singlet-triplet splitting of the former anionic 
groundstate of about 18 meV, and thus a triplet as an¬ 
ionic ground state. 

In this work we add the SOI to our analysis. We find 
that it further removes the triplet degeneracy by inducing 
splittings of few tenths of meV. Moreover, in combina¬ 
tion with exchange correlations, it produces a magnetic 
anisotropy which can in turn be captured by an effective 
spin Hamiltonian. 

In general, the accurate calculation of the many-body 
properties of metallorganic molecules, like the molecular 
magnets or our CuPc, is a highly nontrivial task. In fact, 
the number of their atomic constituents is large enough 
that exact diagonalization is not possible and standard 
density-functional schemes have difficulties in capturing 


short ranged electron-electron correlation^^. In order to 
reduce the size of the many-body Fock space, we use a 
basis of frontier molecular orb itals as the starting point to 
include electronic correlationJ^^^ and construct a gener¬ 
alized Hubbard Hamiltonian. Furthermore, the symme¬ 
try of the molecule greatly helps to reduce the number 
of matrix elements one has to calculate in this basis. 


To probe both SOI induced splittings and magnetic 
anisotropy, we further investigated the current character¬ 
istics of a CuPc molecule in an S TM configuration similar 
to the experiments in Refs.^^^^ the molecule is put on a 
thin insulating layer grown on top of a conducting sub¬ 
strate. The layer functions as a tunneling barrier and 
decouples the molecule from the substrate. Hence the 
CuPc molecule acts as a molecular quantum dot weakly 
coupled by tunneling barriers to metallic leads (here the 
STM tip and the substrate). This quantum dot configu¬ 
ration should be favourable to experimentally probe SOI 
splittings and magnetic anisotropies when an external 
magnetic field is applied to the system, in analogy to the 
experiments in Reffl Indeed, we demonstrate that ex¬ 
perimentally resolvable SOI splitting should be observed 
at magnetic fields of a few Tesla. 


The paper is organized as follows: In Sec. [H] we derive 
a microscopic Hamiltonian for CuPc in the frontier or¬ 
bital basis which includes exchange correlations and the 
SOI. This Hamiltonian is diagonalized exactly and used 
in further spectral analysis and transport calculations. 
Its spectrum is also used to benchmark the prediction 
of an effective spin Hamiltonian which well captures the 
low energy properties of CuPc both in its neutral and an¬ 
ionic configurations. Finally, transport calculations with 
and without magnetic fields are presented and SOI in¬ 
duced signatures are analyzed. Section m contains our 
conclusions. 
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FIG. 1. (a) Geometry and atomic composition of GuPc. (b) 

Single particle energies of relevant molecular orbitals. Black 
(grey) circles depict the tt (ct) character of the corresponding 
orbital. The color (diameter) of the inner circles characterizes 
the type (weight) of the metal orbital contribution on the 
corresponding molecular orbital, (c) Depiction of the four 
frontier orbitals retained in this work: SOMO (S'), HOMO 
{H) and LUMO^^^/y^ {^^zxiyz)- 


II. RESULTS AND DISCUSSION 


A. Microscopic model Hamiltonian for CuPc 


The focus of this section is the establishment of a min¬ 
imal model Hamiltonian for an isolated CuPc molecule 
capable to account for both electron-electron interaction 
and spin-orbit coupling effects. As discussed below, pa¬ 
rameters are fixed such that experimental observations 
for the singlet-triplet splitting as well as positions of 
anionic and cationic resonanceJ^^ are satisfactorily re¬ 
produced. In its most general form and for a generic 
molecule such Hamiltonian reads 


Hmol — Ho + Vee + Vs05 (1) 


where the single-particle Hamiltonian of the molecule is 
given by Hq, Vge describes electronic interactions and 
Vso accounts for the spin-orbit interaction (SOI). 


1. Single particle Hamiltonian for CuPc 

The one-body Hamiltonian Hq, written in the atomic 
basis I a), reads 

Ho = X! + bcp) (2) 

a^a 

where a is a multi-index combining atomic species and 
orbital quantum number at position Tq,, see Fig. 0(a). 
For the ligand we consider the set of all 2s (Is for hydro¬ 
gen), 2px 2p^ orbitals as the cr-system, and conse¬ 
quently the set of 2pz orbitals as the 7r-system. On the 
metal, the Sd^,^, 3 da, 2 _^ 2 , 3(1^2 and 4s orbitals contribute 
to the cr-system, while the Sd;^^; and 3dyz belong to the tt- 
system. This basis yields a total of 195 valence electrons 
for neutral CuPc. Atomic onsite energ ies C g and geomet¬ 
rical parameters were taken from Refs.^^^. The hopping 
matrix elements bap in Eq. 0are obtained by using the 
Slater-Kostei^ and HarrisorP^ LCAO schemes, similar 
to Ref.l^. Numerical diagonalization of Hq finally yields 
single particle energies see Fig. [T|(b ), and molecular 
orbitals |icr) = Cia cf. AppTTO 

Stemming from Hartree-Fock calculations for isolated 
atom^^, the atomic onsite energies do not take into 
account the ionic background of the molecule and crystal 
field contributions. Therefore, molecular orbital energies 
Ci have to be renormalized with parameters Si to coun¬ 
teract this shortage, yielding (cf. App. 

Ho = + (3) 

ia 

In this work we use a constant shift Si = S = l.b eV. 

Due to the odd number of valence electrons, in its neu¬ 
tral configuration CuPc has a singly occupied molecular 
orbital (SOMO). The latter also does not become dou¬ 
bly occupied when the molecule is in its anionic ground- 
stat^. Hence, the orbitals most relevant for transport 
(frontier orbitals) are the SOMO (S'), the HOMO (H) 
and the two degenerate LUMOs {Lzxjyz)', which trans¬ 
form according to the aiu and eg irreducible repre¬ 
sentations of the point group of CuPc ( 04 /^), respectively. 
They are depicted in Fig. (c). The LUMO orbitals in 
their real-valued representations, \Lzx) and \Lyz), have 
equal contributions cp ^ 0.097 on either Sd;^^, and Sd^;^ 
orbitals on the metal, respectively. Due to their degen¬ 
eracy, they can be transformed into their complex, rota¬ 
tional invariant representations: 

\L±)=T2-^^^{\Lzx)±i\Lyz)) 

= T2 / (^\Lzx)pc ^ 

T2 / CL(^\3dzx)cu^'^\^^yz)cu^ 

~ l^^)pc + |3,2, ±l)cu 5 (4) 

where |3,2,±1 )q^ is the n = 3 metal orbital with an¬ 
gular momentum £ = 2 and magnetic quantum number 
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m = zbl. To distinguish contributions from the pure ph- 
thalocyanine (Pc) ligand and the copper (Cu) center, we 
introduced |•)p^ and respectively. Likewise, with 

cs ~ 0.90, we can write for the SOMO: 

1-^) = \/l - 4 |S')p^ +C5 |3da;2_j^2)^^ 

“ \/bT| |S')p^ + 2 C 5 ( |3,2, —2)qjj + |3,2, 2)^^^), 

(5) 

where |3, 2, ±2)^^ is the n = 3 metal orbital with angu¬ 
lar momentum i = 2 and projection m = ±2 onto the 
^-axis. Finally, the HOMO has no metal contributions 
and thus we have trivially \H) = \H)p^. The representa¬ 
tions introduced in Eq. 0 have the advantage that the 
four frontier orbitals can then be characterized by the 
phases (pi acquired under rotations of ^ around the main 
molecular symmetry axis; for the SOMO (ps = tt, for the 


HOMO ph = 0 and for the two LUMOS pl± = • 

2. Many-body Hamiltonian in the frontier orbitals basis 

In order to set up a minimal many-body Hamilto¬ 
nian, we restrict the full Fock space to many-body states 
spanned by the SOMO (S'), the HOMO (H) and the 
two LUMO (T±) orbitals and write Eq. 0 in this basis. 
Hence, for neutral CuPc the number of electrons popu¬ 
lating the frontier orbitals is A^o = 3. 

We exploit the distinct phases acquired by the frontier 
orbitals under 90 degree rotations to determine selection 
rules for the matrix elements Vijki in Vee, 

V- = E E dLdL'd,,,d,.,, (6) 

ijkl era' 

namely Vi^ki 7 ^ 0 if -h 0 /c — = 27r • n, n G Z, cf. 

App. E Equation 0 in this basis then reads 


V 


ee 




\ E 

[ij] 


^ E E H 

[ij] 


^ E E ^ E E (Afe + h.c.) 

[ij] a [ijk] cr 


(7) 


where the indices now run over the set of frontier 

orbitals, and the notation [ijkl] means that the sum runs 
only over unlike indices, i.e. k^ and I are different 

from each other in the corresponding sum. The abbrevi¬ 
ations we introduced in Eq. 0 are the orbital Coulomb 
interaction Ui = Van, the inter-orbital Coulomb interac¬ 
tion Uij = Viijj^ the exchange integral = Vijji^ the 
ordinary pair hopping term Jf- = Vijij and the split pair 
hopping term = Vijik- Contributions with four dif¬ 
ferent indices are found to be very small (on the order 
of yueV) and thus omitted in this work. The matrix ele¬ 
ments Vijki are calculated numerically using Monte Carlo 
integratiorP^ and renormalized with a dielectric constant 
Sr = 2.2 in order to account for screening by frozen or¬ 
bital A table (cf. Tab. 0 with the numerically evalu¬ 
ated interaction constants is found in App 


3. Spin-orbit interaetion (SOI) in the frontier orbitals basis 

A perturbative contribution to the bare one-body 
Hamiltonian Hq relevant in molecular systems is provided 
by the SOI. In the following we derive an effective spin- 
orbit coupling operator acting on the subset of frontier 


orbitals. The atomic SOI operator reads 

VsO = ^ • Sa, (8) 

a,£oc 

where a and run over all atoms and shells, respectively. 
By evaluating Eq. ^ only on the central copper atom, 
i.e. i = 2 and a = Cu, Vso in second quantization is 
given by 

^SO Ccu I ^ ^ 

\m=-2 ^ 

+ + di^dg^ + h.C.^ 

+ + h.c.^ ^ , (9) 

where dj^^ creates an electron with spin a on the cop¬ 
per atom in the orbital specified by {£ = 2,m). For an 
electron in the 3d-shell of Cu we use ^cu ~ 100 me\^. 
Projecting Eq. 0 onto the minimal set of frontier or¬ 
bitals then yields: 

Vso =Ai ^ r — d^ri^Lri) 
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FIG. 2. Lowest lying anionic states of CuPc, together with 
their grade of degeneracy d. Without exchange and SOI, the 
anionic groundstate is eightfold degenerate. When exchange 
interaction between SOMO and LUMOs is introduced, the 
degeneracy is lifted, yielding two triplets and two singlets be¬ 
cause of the orbital degeneracy of the LUMO. SOI further 
splits the triplet states, generating a twofold degenerate an¬ 
ionic groundstate consisting of the states and T^. 

+A 2 -h + h-c.^ , (10) 

where Ai = |clP = 0.47 meV and A 2 = Ccu = 
6.16 meV are now effective spin-orbit coupling constants. 
A similar analysis of SOI in CuPc, laying more focus on 
the central Cu atom, can be found in Ref.^^. 

Finally, many body eigenenergies E^k and eigenstates 
\Nk), labelled after particle number N and state index 
k, are obtained by exact numerical diagonalization of 
Hmoi in the frontier orbitals basis. Despite numerically 
tractable, the problem described by Hmoi is still highly 
intricate, as the Fock space has dimension 4^ = 256. In 
reality, though, only few low-lying many-body states are 
relevant at low energies, what enables further simplifi¬ 
cation and even an analytical treatment, as discussed in 
the next subsection. 


B. Low-energy spectrum of CuPc and effective 
spin Hamiltonian 

In the following we will analyze the neutral and anionic 
low-energy part of the many-body spectrum of CuPc and 
establish an effective Hamiltonian which enables us to 
analyze the low-energy behaviour in a more lucid way. 
To this extent, we start by observing that Hmoi (in the 
considered particle number subblocks) contains different 
energy scales, in particular, U > J > A, what suggests 
a hierarchy of steps. We use U, J and A to denote the 
set of all Hubbard-like parameters (f/^, Uij)^ all exchange 
parameters {JifiJfjiJ'ijk) nnd all SOI parameters (A^), 
respectively. As a first step we set both the exchange (J) 
and SOI (A) contributions to Hmoi to zero and determine 
the neutral and anionic groundstates. In a second and 
third step exchange and SOI are added, respectively. 


1. Neutral low-energy speetrum 

In the neutral low-energy part of the spectrum, we 
retain the two spin-degenerate groundstates of HmoK'/ = 
0,A = 0), 

\No,a):=dljn), ( 11 ) 

with corresponding energy Ef^^. Here we defined \ ft) = 
|0). The groundstates in Eq. ( pTj ) are neither af¬ 
fected by Vso nor by the exchange terms in Eq. Q . Triv¬ 
ially, the effective Hamiltonian in the basis of\No^ga) 
reads: 

( 12 ) 

In principle Eq. 0 also contains terms which act on 
the neutral groundstate, like for example pair hopping 
terms proportional to cause admixtures 

with other many-body states. However, according to our 
full numerical calculations, these admixtures are rather 
small and do not affect transitions between neutral and 
anionic states. 


2. Anionie low-energy speetrum 

Continuing with the anionic low-energy part of the 
spectrum of Hmoi(^ = 0, A = 0), we find an eightfold 
degenerate groundstate: 

\No + 1, TCTcr') := dkdba' \^) > (13) 

with corresponding energy The eightfold degen¬ 

eracy comes from the two unpaired spins in either SOMO 
or LUMO and the orbital degeneracy of the LUMO or¬ 
bitals. In order to make the anionic eigenstates also 
eigenstates of the spin operators and Sz^ they can 
be rewritten as 

l^r) = 1 ^) ’ 

|Tt) = d^t^Lt 111). 

|Tt) = (dstdiri + ds^diTt) |11) ’ 

|T;)=d^^dU|^^). (14) 

The orbital degeneracy of the LUMOs, expressed by the 
index r, is responsible for the two sets of singlets (total 
spin S = 0) and triplets (total spin S = 1). Considering 
exchange interaction in a second step, we find that only 
the term in Eq. Q, 

“ ^ ^Tl (^^Sa^Lra “ 5 ( 15 ) 

rcr 

directly determines the low-energy structure of the an¬ 
ionic low-energy part because of the singly occupied 
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SOMO and LUMOs: The degeneracy between singlets 
and triplets is lifted, see Fig. and we obtain 

Es = 

Et = E%^ - (16) 

for the singlets and triplets, respectively. 

Finally, to analyze in a third step how Vso affects the 
low-energy part of the anionic part of the spectrum, in 
particular which degeneracies are lifted, we treat it as a 
perturbation and apply second order perturbation theory 
to obtain the energy shifts. To this end, some additional 
states have to be considerd. They are listed in App. 0 
The states and experience a downshift due to 
Vso and become the groundstates. Measuring energies 
with respect to we get 

2 \2 

AE = AE^^ = -Ai - (17) 

+ - -h 

see Fig. Note that in our numerical calculations 
and are mixed and the degeneracy of the resulting 

states is lifted by a small shift in the range of some /ieV. 

A more detailed discussion concerning the mixing of 
and can be found in App. The next states are 
and with 

\2 \2 

AErj^O = AEr^O = -^ ?-. (18) 

+ - 2J- Ai-J- ^ ^ 

Due to their quadratic dependence on Ai and A 2 , these 
states change very little with Vso- The degeneracy of 
the states and TI is lifted by the mixing of these 
states through Vso- We find 

l<»> = ^(|3’^) + |3’r)), (19) 

l« = ^(|T)()-|Tr)), (20) 



B.(T) 



B.(T) 


FIG. 3. (a) Dependence of the single particle orbital energies 
on the magnetic field strength. From this, the effective or¬ 
bital moment of the LUMOs, here depicted in their complex 
representation, can be extracted as /Xorb = 33.7 /xeVT”^. The 
energies of the SOMO and HOMO orbitals depend quadrat- 
ically on the magnetic field and involve a much lower scale 
than the LUMOs, as seen in the close-up in panel (b). 


Equation (24) is one major result of this work. It shows 
that, similar to the well studied molecular magnet^^M^, 
the interplay of spin-orbit coupling and exchange inter¬ 
actions yield magnetic anisotropies which can be cap¬ 
tured by effective spin Hamiltonians. Noticeably, be¬ 
cause Eq. (24) was derived from the microscopic molec¬ 
ular Hamiltonian HmoU if was possible to check that de¬ 
viations are in the /ieV range and only of quantitative 
nature by comparison of the spectrum to the numerically 
evaluated one. 


C. Interaction with magnetic fields 

An experimentally accessible way to probe magnetic 
anisotropies is to apply external magnetic fields. In or¬ 
der to account for interactions of orbitals with magnetic 
fields, the atomic hopping matrix elements in Eq. (© 
have to be corrected with Peierls phase factors. 


where for \P) we omitted smaller additional contributions 
from other states. The energies change according to 


AE{a) = X,, 


Tex 

^SL 


( 21 ) 

( 22 ) 


Eor further details we refer to App. Einally, the sin¬ 
glets S+ and S_, similar to and T^, change very 
little (with respect to Es)' 


AEs, 


Af _ Aj 

2'^Il — 'A|l 


(23) 


By introducing f := an approximate Hamil¬ 

tonian up to first order in Vso can be given for the Vq + I 
particle subblock: 


- JTl (52 - 1) + Ai tS, . (24) 


(25) 

where, using the gauge A = —Bzyx, the phase is given 
by 


0a/3 — 2^ T y^) {Xcx xp) . (26) 

Here {x^^ya) are the in-plane atomic coordinates. Ow¬ 
ing to the planar geometry of CuPc, depends only on 
the ^-component Bz of the magnetic field B. In Eig. 
we show the dependence of the energies of the frontier 
molecular orbitals on the strength of the magnetic field 
in ^(-direction, Bz. Eor the two LUMOs we observe a lin¬ 
ear dependence on the magnetic field, yielding an effec¬ 
tive orbital moment of /iorb = 33.7 yeVT~^. Hereby the 
LUMO—(+) goes down (up) in energy with Bz^ see Eig.[^ 
(a). The energies of the HOMO and the SOMO however 
scale quadratically with the magnetic field at a much 
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lower scale, cf. Fig. (b). This behaviour is expected, 
since the aiu and big representations have characters +1 
under rotations, which transform to Thus 

the energies of HOMO and SOMO can not depend on the 
sign of Bz and must move at least quadratically with B^. 
The two-dimensional eg representation on the other hand 
has zero character under C 2 rotations, which implies that 
the constituents of eg transform under such rotations ei¬ 
ther with different signs or into each other; indeed under 
a C 2 rotation LUMO+ is mapped onto LUMO— and vice 
versa. 

Finally, the interaction of electronic spins with mag¬ 
netic fields is represented by adding a Zeeman term Vz 
to Eq. Q, 


Hmol ^ Hmol + Vz = Hmol + 5 'S'Mb S • B, (27) 

where gs = 2 and S is the total spin operator on the 
molecule written in the frontier orbital basis. 


1. Effective low-energy Hamiltonian 


Putting everything together, an effective low-energy 
Hamiltonian including magnetic interaction terms for 
both orbital and spin degrees of freedom can thus be 
given. It reads 

+fl5MBS ■ B, (28) 


N 


where H^ 
low-energy 
and (24). 


is the Hamiltonian for the corresponding 
V-particle subblock as given by Eqs. (12) 


D. Dynamics and transport 

1. Reduced density operator and current 

The transport calculations for the molecule in an STM 
setup are done by using the formalism introduced in ear¬ 
lier work4^ ^^ l ^^ ( Eor the sake of clarity, in the following 
we briefly discuss the main steps to obtain the current 
through the molecule. The full system is described by 
the Hamiltonian 


H — Hmol + Hie + Hs + Ht + Htun, (29) 

where Hmol describes the isolated molecule, see Eq. §■ 
To incorporate image charge effects in our model, leading 
to renormalizations of the energies of the systems charged 
stated, we included a term HiJ^, 

Hic = -(5ic(iV-7Vo)', (30) 

where N is the particle number operator on the molecule. 
Electrostatic considerations regarding the geometry of 



1 ^ 
¥ 


-1 


14 (V) 


FIG. 4. Current and differential conductance curves exhibit¬ 
ing the anionic (cationic) resonance at positive (negative) bias 
voltage. Note that in contrast to all other results in this work, 
this curve is taken at a temperature of 60 K to emphasize the 
resonances in the dlf dV curve. 


the STM setup yielded Sic ~ 0-3 e\^. The Hamilto¬ 
nians Hs and Ht corresponding to substrate (S) and tip 
(T), respectively, are describing noninteracting electronic 
leads. They read 

Hr7=S,T = ^ 

kcr 

where creates an electron in lead g with spin a and 

momentum k. The tunneling Hamiltonian Htun finally is 
given by 

Htun = cj^kcr^icr + h.C.. (32) 

?7kicr 

It contains the tunneling matrix elements t^-, which are 
obtained by calculating the overlap between the lead 
wavefunctions \gk.) and the molecular orbitals 

Finally, the dynamics of the transport itself is calcu¬ 
lated by evaluating the generalized master equation, 

Pred — ^^[Pred]^ (33) 

for the reduced density operatoJ^^^ pj-ed = Trs,T(p)- 
The Liouvillian superoperator 

£ = £s + £t + £rel (34) 

contains the terms Cs and £t describing tunneling from 
and to the substrate and the tip, respectively. To ac¬ 
count for relaxation processes leading to de-excitation of 
molecular excited states, we included a relaxation term 
£reb analogously to Ref.^ 

[p] = -l(^p- ptE \Nk) {Nk \E ■ (35) 
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It depends on the deviation of p from the thermal solution 
which is given by a Boltzmann distribution: 


P 


th 


E 


e-P^Ni 


\Nk) {Nk \, 


(36) 


with [3 = (/cbT)~^. Since C^ei acts separately on each 
A/'-particle subblock, it conserves the particle number on 
the molecule and thus does not contribute to transport 
directly. In this work, the relaxation factor ^ is around 
the same order of magnitude as the mean tip tunneling 
rate onto the molecule. In particular, we are interested in 
the stationary solution for which p^^ = E[p^^] = 0. 
Finally, the current through the system in the stationary 
limit can be evaluated as 

(Is + ^t) = ^ (iV) = IV^oi {Ncm) = 0, (37) 

yielding the current operator for lead p as /^ = 


2. Transport characteristics 


In this work, a tip-molecule distance of 5 A was used 
and simulations were done at the temperature T = 1 K. 
We assumed a substrate workfunction of 0o = 4.65 eV 
and a renormalization of the single particle energies 
Si = 6 = 1.5 eV (cf. Eq. Q). Numerical results for 
the current and the differential conductance, according to 
Eq. (37) and using the full Hamiltonian Hmoi in Eq- (29), 
are shown in Eig. (|^. Anionic (cationic) resonances at 
positive (negative) bias voltages are clearly seen. We find 
a very good agreement between our numerically evalu¬ 
ated positions of the cationic and anionic resonances with 
those of the experiment in Ref.t^^, where a Cu(lOO) sub¬ 
strate was used. 

Notice that, in our model, the bias voltage at which 
a tip-mediated transition from the mth neutral state to 
the nth anionic state of the molecule is happening is 


Hr. 


1 


aT\e\ 


(-Ewo + ljn ^No,m ^ic “ 1 “ 0o) 5 (^^) 


where e is the electron charge and accounts for the 
fact that in STM setups the bias voltage drops asymetri- 
cally across the junction. We are usin g a x = 0.59 for the 
tip and as = —0.16 for the substrate^. If given with¬ 
out indices, Ides denotes the bias voltage corresponding 
to the groundstate-to-groundstate resonance. 

The negative differential conductance at large negative 
bias in Eig. is caused by blocking due to population of 
excited states of the molecule. This has already been dis¬ 
cussed in some previous worlP^ and will not be of further 
interest here. 


3. Transport simulations at finite magnetic fields 

In Eig.j^we show the splitting of the anionic resonance 
with applied magnetic field in a dl/dV map. In the up¬ 



B(T) 


FIG. 5. Differential conductance maps as a function of the 
strength Bz of the magnetic field in z-direction. Upper (lower) 
panel: Spin-orbit interaction switched off (on). Solid and 
dashed lines depict the addition spect rum as calculated from 
the effective spin Hamiltonian, cf. Eq. ( |28| ). Transitions start¬ 
ing from the neutral groundstate are denoted by solid lines, 
those from the neutral excited state by dashed lines. 


per panel SOI is switched off, whereas in the lower panel 
it is switched on. One striking difference at first glance 
is the zero-field splitting for nonvanishing SOI, which is 
proportional to Ai but enhanced by the bias drop, cf. 
Eq. (38). Eor vanishing SOI, when Sz is a good quan¬ 
tum number, we can readily identify the corresponding 
transitions by using the effective spin Hamiltonian intro¬ 
duced in Eq. (28). In the following, transitions from the 
neutral groundstate will be denoted by arable numbers: 


(1) ; |iVo4)^|T:) 

(2) : |iVo,;)^|T;) 

(3) ; |iVo,;)^|T°) 

(4) : |iVo,;)^in), 

while transitions from the neutral excited state will be 
denoted by Roman numerals: 

(i): |Aro,t)^|T°) 
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FIG. 6. Differential conductance maps vs. the angle 0, formed by the applied magnetic field with the z-axis. Left (right) 
panels are without (with) SOI. Upper, middle and lower panels are calculated for a magnetic held strength of 1 T, 3 T and 
8 T, respectively. Solid and dashed lines depict the addition spectrum as calculated from the effective spin Hamiltonian, cf. 
Eq. (28). Transitions starting from the neutral groundstate are denoted by solid lines, those from the neutral excited state by 
dashed lines. 


(ii) ; 

|iVo,t)^ 

|T°+) 

(Hi) : 

\NoA)^ 

|Tl) 

(iv) : 

|iVo,t)^ 

lU) 


Other transitions are forbidden due to the selection rule 
for Sz, ASz = The reason for the splitting into four 
lines observed in Fig. top is that the orbital moment of 


the LUMO is not of the same size as the Bohr magneton. 

For nonvanishing SOI, see lower panel of Fig. the 
definite assignment of transitions is not straightforward, 
at least for small magnetic fields. Since and are 
shifted downward by SOI, transition (2) now is the lowest 
lying transition, whereas transition (1) is shifted upward 
due to the positive contribution +Ai to TI. Further- 
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more, transition {iv) is the only excited state transition 
which can be definitely assigned to a line in the lower 
panel in Fig. 

Figure finally shows dl/dV maps as a function of 
the angle 0 between the magnetic field and the z-axis. 
Hereby panels (a), (b) and (c) show results obtained with 
vanishing SOI and panels (d), (e) and (f) are for finite 
SOL Again, the results were fitted using the effective spin 
Hamiltonian introduced in Eq. ( [2^ with good agreement. 
The respective transitions can be identified by checking 
the assigned transitions in Fig. at the corresponding 
field strength. 

Already at |B| = H = 1 T, cf. (a) and (d), the in¬ 
fluence of SOI can be clearly seen. While for vanishing 
SOI any anisotropy of the dl/dV map is hidden beneath 
the temperature broadening, for finite SOI a slight 0- 
dependence can be observed. For H = 3 T, now also in 
the vanishing SOI case, Fig.|^(b), a slight anisotropy due 
to the orbital moment of the LUMOs can be observed, 
although still blurred by temperature. Again, at finite 
SOI in Fig. (e) there is a much more pronounced de¬ 
pendence on 0. The high conductance areas at ^ = 0° 
and 0 = 180° for W — Kes ~ 0.8 meV correspond to 
the high conductance area in the middle of Fig. [^bot¬ 
tom, where many transitions are taking place at the same 
time. At H = 8 T, the magnetic field is dominating and 
a characteristic double cosine-like behaviour of the reso¬ 
nances can be observed, for both the case with no SOI, 
Fig. [^(c), and finite SOI, Fig. |^(f). For vanishing SOI, 
this behaviour is caused by the orbital moment of the 
LUMOs, since they interchange their positions when go¬ 
ing from Bz to —Bz. The overall splitting between the 
double cosines, most evident at 6> = 90°, is caused by 
the Zeeman term. The results for H = 8 T in Fig. 

(f) at finite SOI are similar to those in Fig. (c), with 
the only difference that the cosine at large biases is more 
stretched, the one at low bias more compressed. 


III. CONCLUSIONS 

We established a model Hamiltonian for CuPc which 
accounts for electron-electron, spin-orbit and magnetic 
interactions in a minimal single particle basis represented 
by four frontier orbitals; the SOMO, the HOMO and 
two degenerate LUMOs. The distinct properties of these 
orbitals under rotations allowed us to deduce selection 
rules for matrix elements of the Coulomb interaction, 
which drastically reduce the number of nonvanishing 
terms and simplify the numerical diagonalization of the 
full many-body Hamiltonian. For the low-energy parts 
of the neutral and anionic blocks of the many-body spec¬ 
trum we could further derive an effective spin Hamilto¬ 
nian, capturing both SOI induced splittings and magnetic 
anisotropy. In order to study fingerprints of the SOI un¬ 
der realistic experimental conditions, we have studied the 
magnetotransport characteristics of a CuPc based junc¬ 
tion in an STM setup. To this extent, a generalized mas¬ 


ter equation for the reduced density matrix associated to 
the full many-body Hamiltonian had to be solved in order 
to numerically obtain both the current and the differen¬ 
tial conductance. Noticeably, by using the effective spin 
Hamiltonian, it was possible to reconstruct the nature 
of the many-body resonances observed in the numerical 
calculations. 

In summary, we believe that our work significantly ad¬ 
vances the present understanding of spin properties of 
CuPc; moreover, the flexibility of our model Hamiltonian 
approach opens new perspectives for the investigation of 
other configurationally similar metallorganic compounds. 
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Appendix A: Transformation from the atomic to the 
molecular orbital basis 


The Hamiltonian of a molecule in Born-Oppenheimer 
approximation, after dismissing terms which only depend 
on the positions of the nuclei and are therefore constant, 
can be written as 


ajSa 

mn 



E 


t y ^ ^ 


aa' 

mnop 


(Al) 


where creates an electron in the atomic orbital 

\ama) with orbital quantum number m and spin a cen¬ 
tered at atom a. Further we have defined 


ham,^n •— T ^am,l3n’! 

where Cam is the energy of orbital m on atom a and 
ham,^n is the hopping integral between orbital m on atom 
a and orbital n on atom /3. All non-hopping terms can 
be condensed in the crystal field correction 

W^m,/ 3 n := E {ama\\^\(ina), (A3) 

7 

wher e is the atomic core potential at r^. Equa¬ 
tion ( |A3[ ) defines the crystal field correction to the single 
particle Hamiltonian. Einally, we have the ordinary ma¬ 
trix elements of the Coulomb interaction. 

The ham,i3n are elements of a matrix h which corre¬ 
sponds to the single particle Hamiltonian of the molecule 
with only onsite energies and hopping terms. After per¬ 
forming a transformation to the molecular orbital basis. 
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in which h is diagonal, |icr) = ^iam \(^ma)^ and us¬ 
ing the approximation that the basis \ama) is orthogonal, 
the Hamiltonian reads: 

ija 

^ijkl ^ia^ka'^la'^ja^ 

ijkl <j(j' 

where now is 

•'j 

= E (A5) 

al3 


One huge simplification which is possible in the molec¬ 
ular orbital basis is the reduction of the size of our Hilbert 
space H, which occurs by retaining few relevant molecu¬ 
lar orbitals only. To this end we split the full molecular 
basis into frozen and dynamic orbitals, where Nf of the 
frozen orbitals are assumed to be always fully occupied 
and the remaining set to be always empty. We do not 
make any assumption about the occupation of the 
dynamic states. Whether these frontier orbitals are 
full or empty depends on the electrochemical potential of 
the molecule, and on whether an exchange of electrons 
with the environment is possible. 

In the occupation number representation a general 
state of the Fock space then looks like 


Appendix B: Symmetries in the frontier orbitals 
basis 


|T) |11... 11) (g) luk^Uki ... ni^nii) (g) |00 ... 00). (Bl) 


2Nf 


2Nd 


2Ne 


Us 

11.352 eV 

Tex _ tP 

^HL — -^HL + L- 

548 meV 

Uh 

1.752 eV 

rex 

^L + L- 

258 meV 

Ul = Ul+l- 

1.808 eV 

^L + L- 

168 meV 

USH 

1.777 eV 

II 

1 

+ 

1 

9 meV 

UsL 

1.993 eV 

II 

2 meV 

Uhl 

1.758 eV 




TABLE I. Major nonvanishing Coulomb integrals between the 
SOMO(5'), the HOMO(H), the LUMO+ and the LUMO-. 
When the LUMOs need to be distinguished, they are denoted 
as L+ or L—, otherwise just by L. All values are calculated 
numerically using Monte Carlo integratiorP^ of the real space 
orbitals depicted in Figs.[^c) and[^ respectively, and renor¬ 
malized by a constant Cr = 2.2. 


In this work we assume the molecule to be neutral un¬ 
der equilibrium conditions, with 195 valence electrons. 
Thus, the orbitals we choose to build up the subspace of 
dynamic orbitals are orbitals Nrs. 97-100, see Fig. [ 11 (b). 
This choice results in the lowest 96 molecular orbitals 
being doubly filled. Note that the choice of the LUMO 
states I/± rather than L^x/yz is convenient due to the 
fact that these orbitals acquire a definite phase upon ro¬ 
tations of 90 degrees around the main symmetry axis of 
the molecule. Specifically, for the four frontier orbitals S', 
H and LT, the acquired phases are (ps = tt, <Ph = 0 and 
(j)L± = ±f, respectively. This in turn imposes symmetry 
constraints on the Hamiltonian (A4). Consider e.g. the 
Coulomb interaction 


yijki = Jj dVi dV2t/>*(ri)V;j(ri)|—L_^*(r2)t/,,(r2). 


(B2) 


Then, in the frontier orbital basis it holds that: 

Vijki = Vijki- (B3) 

Therefore a given matrix element of the Coulomb in¬ 
teraction Vijki is different from zero only if the sum of 
the corresponding phases adds up to multiples of 27r: 
4>i — + 0/e — = 27r • n, n G Z. In Tab. we list 

all nonvanishing matrix elements of the Coulomb inter¬ 
action which are used in this work. For the crystal field 
correction h can be shown that: 

f'j 

AF/y = AV/y (B4) 

^AlAf = Al/A°"%, (B5) 

since all phases 0^ are different; 0^ 7 ^ 0^ for i ^ j. Hence 


AH/2^ is diagonal in the basis. In the follow¬ 

ing we treat the AH/°^ as free parameters and include 
them in the paramteter 5i entering Eq. (§. 


Appendix C: Details on the perturbative treatment 
of SOI 


In addition to the states introduced in Eq. (14), the 
following states must be also taken into account when 
performing second order perturbation theory: 


\Lt t,LTi) = dVdV\n), 
\Lt(J, Lfaj = dfradLa' 1^) > 
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(Cl) 


with = ELTa,Lfa' = and Es^,si = A 2 . In 

the basis introduced in Eqs. (14) and (iCll), Vso is block- 


diagonal and decomposes into six subblocks: two three- 
dimensional, two two-dimensional, one four-dimensional 
and one one-dimensional subblocks. 

The four dimensional subblock describes the effects 
of SOI on the T| and TI states. Written in the ba¬ 
sis {|T^),|TI),|I/+ t: S i)}: Hamiltonian 

reads 


/-m 0 0 0 \ 

I 0 - 0 0 I 

0 0 Ai 0 

\ 0 0 0 A 2 / 

/ Ai 0 -V 2 A 2 

0 Ai V2X2 

—\pl\2 '/2X2 Ai 

V V2X2 -V2X2 0 

The degeneracy of the unperturbed states and TI 
and the fact that there are no matrix-elements which cou¬ 
ple these states require the use of second order degener¬ 
ate perturbation theory. Applying it yields the following 
matrix M: 


V 2 X 2 \ 

-^^2 _ (C 2 ) 

0 / 


M = A- 


1 -1 

-1 1 


where the prefactor A is given by 

A = -2A2' ^ 


+ ^SL ^2 + JsL 


(C3) 


(C4) 


Diagonalization of M gives the second-order energy cor¬ 
rections 


AE{a) = Ai, 
AE{I3) = Ai - 4A2 


Ai + A 2 + Jg\ 


(C5) 
, (C6) 


and the correct linear combinations of the states T^ and 
TI: 


H = ^(|Tp + |TI)) (C7) 

l/3) = ^(|Tp-|TI)). (C8) 

Writing H in the basis {|(a),|;d),|I/+5'|)} 
yields: 


/-JFl 

0 

0 

V 0 


0 0 
-JTl 0 
0 Ai 
0 0 


0 

0 

0 

A 2 


/Ai 0 0 0 \ 

0 Ai —2A2 2A2 

0 — 2 A 2 Ai 0 

Vo 2 A 2 0 0 / 


(C9) 


We see that ja) stays unaffected by the perturbation, 
whereas \P) will change: 


1/3) ^1/3 )+2 \L+ t,L- I) 

- 2 —A—IS't, 5 ;). (CIO) 

^2 + 

The mixing of and is caused by a pair-hopping 
term in the Hamiltonian, more precisely by 

\Jl+L- E {^L+AL+iT^L-iT^L-,T + h-C.) , (Cll) 

a 


which couples and to the following states: 

1 ^) = 1 ^) ’ 

1^) = (dL+t^L+^ + |0) 5 (C12) 

with corresponding energies Ea and = Ea + 

Then, after introducing 

|Ti) = T(|t;) + |tP), 

IT 2 ) = T (|t;) - |Tt)), (C13) 

the Hamiltonian in the basis of these four states can be 
written as 

" = ( 1 ;‘ h 1 ) • (cn) 

with 

"““(■■'“aI"' £() 

and 

=(“■'"aI k) ^ 

Diagonalization finally yields the four states 

|1) = ^^(|Ti)+76|6)), 

Vl-7b 

|2) = ^A=(|T2)+7a|a)), 

Vl-7a 

|i) = ^A==(|6)-7b|Ti)), 

V 1 - % 
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|2)= ^-^ (|a)-7a|T2)), (C17) 

Vl-7a 

with the admixture 7 „/(| « f • Their energies are 

' ‘^a/h-r-JsL 

approximately 


El 

E 2 


Eh 


Ea + 


Ai 


Eh + ^ 

Ea + Js\ E Ai 


(C18) 


El « 

1 

1 

Ai 


Eb + 

+ Ai^ 

E2- 

1 

1 

Ai 


Ea + Js\ 

+ Ai 


This analysis reproduces mixing and energy splittings 
consistent with our numerical calculations. 
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